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The stationary points of the potential energy function of the model on a two-dimensional 
square lattice with nearest-neighbor interactions are studied by means of two numerical methods: a 
numerical homotopy continuation method and a globaUy-convergent Newton-Raphson method. We 
analyze the properties of the stationary points, in particular with respect to a number of quantities 
that have been conjectured to display signatures of the thermodynamic phase transition of the model. 
Although no such signatures are found for the nearest-neighbor i^* model, our study illustrates the 
strengths and weaknesses of the numerical methods employed. 
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I. INTRODUCTION 

The stationary points of the potential energy function 
or other classical energy functions can be employed to 
calculate or estimate certain physical quantities of inter- 
est. Well-known examples include transition state theory 
or Kramers's reaction rate theory for the thermally ac- 
tivated escape from metastable states, where the barrier 
height (corresponding to the difference between poten- 
tial energies at certain stationary points of the potential 
energy function) plays an essential role. More recently, 
a large variety of related techniques has become popular 
under the name of energy landscape methods [1], allow- 
ing for applications to many-body systems as diverse as 
metallic clusters, biomolecules and their folding transi- 
tions, or glass formers undergoing a glass transition. 

In the late 1990s it was observed that properties of 
stationary points of the potential energy function V , i. e. 
configuration space points satisfying dV{q^) = 0, re- 
flect in dynamical and statistical physical quantities si- 
multaneously and show pronounced signatures near a 
phase transition [2] . This observation sparked quite some 
research activity, reviewed in [3] , including a theorem by 
Franzosi and Pettini asserting that, at least for a certain 
class of models, stationary points with V{q^)/N = Vc 
are indispensable for the occurrence of an equilibrium 
phase transition at potential energy [4] . This theorem 
requires a number of conditions to be satisfied: The po- 
tential energy function V has to be smooth, confining, 
and of short-range (sec [4] for a complete list of con- 
ditions and their definitions). At the time when these 
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papers were published, one might have still hoped that 
some of the conditions on V were merely technical, but 
not essential for the result. However, it became clear 
soon that the result can not be extended to long-range 
interacting models [5] , nor to non-confining potentials [6] : 
These classes of potentials comprise cases which are par- 
ticularly amenable to analytic calculations, and a direct 
relation between phase transitions and stationary points 
of V could be ruled out through exactly solvable coun- 
terexamples. 

Originally, the incentive for the study reported in the 
present article was to investigate the stationary points 
of a model that satisfies all the conditions required by 
Franzosi and Pettini [4]. This is not an easy task, as in 
this class there are no exactly solvable models which have 
a phase transition [7] . As a model to study, we then opted 
for the nearest-neighbor (j)^ model on a two-dimensional 
square lattice. This model, though not exactly solvable, 
appears to be relatively simple. Moreover, results on the 
stationary points of its long-range version were known 
and readily available for comparison [5]. 

Much to our surprise, we found that all stationary 
points q^ of the potential energy function V have non- 
positive potential energies, i.e., V{q^) < 0. From this 
observation, one can conclude that the result of Franzosi 
and Pettini, allegedly proven in [4], is false. Further- 
more, a numerical method put forward in [8] and applied 
to the very same two-dimensional cj)'^ model yields incor- 
rect results. These findings, and a discussion of their im- 
plications, have been published in a Letter [9]. The non- 
positivity of the stationary energies V{q^) was established 
in that Letter analytically, supported by results obtained 
with two different numerical methods. The main pur- 
pose of the present article is to give a detailed account of 
these numerical methods and to present a more detailed 
analysis of the properties of the stationary points of the 
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two-dimensional nearest-neighbor cj)'^ model. 

In Sec. II, this model is introduced and some of its 
thermodynamic properties are reviewed. In Sec. Ill, the 
first of the numerical methods, namely homotopy contin- 
uation, is discussed. It is an algebraic-geometrical tech- 
nique devised to obtain all isolated stationary points of 
a given system of multivariate polynomial equations, but 
is restricted to fairly small lattice sizes. We have applied 
this method to square lattices of sizes 3x3 and 4x4. 
The stationary points obtained are analyzed with respect 
to their number, potential energies, indices, and Hessian 
determinants in Sec. IV. The second numerical method, 
discussed in Sec. V, makes use of a globally convergent 
version of the Newton- Raphson algorithm for searching 
the zeros of a real-valued function. It can be applied to 
larger lattice sizes, but provides in general only a subset 
of the stationary points. We summarize and discuss our 
findings in the concluding Sec. VI. 



II. TWO-DIMENSIONAL 
NEAREST-NEIGHBOR (j)'* MODEL 

On a finite square lattice A C consisting oi N = 
sites, a real degree of freedom (j>i is assigned to each lattice 
site i € A. By we denote the subset of A consisting 
of the four nearest-neighboring sites of i on the lattice 
under the assumption of periodic boundary conditions. 
The potential energy function of this model is given by 
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where q = (qi, . . . tQn) denotes a point in configuration 
space r = K,^ [10]. The parameter J > determines the 
coupling strength between nearest-neighboring sites and 
the parameters A, /x > characterize a local double-well 
potential each degree of freedom is experiencing. 

In the thermodynamic limit N ^ oo this model is 
known to undergo, at some critical temperature Tc, a 
continuous phase transition, in the sense that the config- 
urational canonical free energy 
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is nonanalytic aX T — T^. The transition is from a "ferro- 
magnetic" phase with nonzero average particle displace- 
ment to a "paramagnetic" phase with vanishing aver- 
age displacement (see [11] for more details as well as for 
Monte Carlo results). 

Since we are interested in whether, and how, the phase 
transition reflects in the properties of the potential en- 
ergy landscape, it is more adequate for our purposes to 
compare not to Tc, but to the critical potential energy 
per lattice site, Vc, of the transition [12]. Both quantities 
are unambiguously related to each other in the thermo- 
dynamic limit via the caloric curve v{T). This is true 



independently of the statistical ensemble used, as these 
ensembles are known to be equivalent for short-range 
models like the one we are studying [13]. 

The critical potential energy Vc is less frequently stud- 
ied, in fact the only data we could find in the literature 
are from Monte Carlo simulations of fairly small system 
sizes = 20 X 20 in [14], with parameter values A = 3/5, 
= 2, and J = 1. We use the same values of A and 
/i^ in the following, but will show results for a range of 
couplings J. Since the value of Vc is a crucial benchmark 
when relating our stationary point analysis to the phase 
transition of the 0^ model, we have performed standard 
Metropolis Monte Carlo simulations for somewhat larger 
system sizes up to 128 x 128 and 10^ lattice sweeps. 

Some of the Monte Carlo results have already been 
reported in [9]. From these plots one can read off a crit- 
ical potential energy per lattice site of roughly Vc ~ 2.2 
for coupling J = 1. An improved estimate could be ob- 
tained by more extensive Monte Carlo simulations and/or 
a finite-size scaling analysis of the data, but the results 
as they are will be sufiicient for our purposes. We have 
determined Vc also for a few other couplings, obtaining 
the rough estimates Vc ~ —5.0 for J = 0.2 and Vc ~ —3.0 
for J = 0.35. 



III. NUMERICAL POLYNOMIAL HOMOTOPY 
CONTINUATION METHOD 

The idea behind numerical continuation methods is to 
first find the solutions of a simple system of equations 
which shares several important features with the given 
system. Then, in a second step, starting from these so- 
lutions one continues them towards the given system in 
a systematic way. Homotopy continuation methods have 
been around already for several decades [15, 16]. With 
more recent machinery like the numerical polynomial ho- 
motopy continuation (NPHC) method used in the present 
article, the method is guaranteed to find all isolated so- 
lutions of systems of polynomial equations [17, 18]. 

We consider a system of m polynomial equations 
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in the variables q = (qi, . . . , qm)'^ , and we assume that all 
solutions of (3) are isolated. Then Bezout's Theorem (see 
Chapter 8 of [17]) asserts that a system of m polynomial 
equations in m variables has at most Oi=i '^i isolated 
solutions where di is the degree of the ith polynomial. 
This bound is called the classical Bezout bound, and it is 
known to be sharp for generic systems [i.e., for generic 
values of the coefficients of the polynomials Pi{q)]. 

The continuation of solutions is formally described by 
the homotopy 



H{q,t)^Piq){l-t)+jtSiq), 
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where 7 is a complex number and 



S{q) 



(5) 



is again a system of m polynomial equations. Varying 
the parameter t G [0,1], H can be deformed from the 
start system H{q, 1) = jS{q) at t ~ 1 into the polynomial 
system of interest, H{q, 0) = P{q) at t = 0. The following 
conditions have to be satisfied in order to guarantee that 
all solutions of P can be computed from this homotopy: 

1. The solutions of S{q) = can be computed. 

2. The number of solutions of S{q) = satisfies the 
classical Bezout bound for P{q) = as an equality. 

3. The solution set of H{q, t) = for t G (0, 1] consists 
of a finite number of smooth paths, called homotopy 
paths, which are parameterized by t. 

4. Every isolated solution of H{q, 0) = P{q) = can 
be reached by some path originating at a solution 
ofH{q,l) = ^S{q) = 0. 

Satisfying the first two criteria hinges on a suitable choice 
of the start system S. Criteria 3 and 4 are guaranteed to 
be satisfied based on the gencricity of the constant 7 in 
(4). Theorem 8.4.1 of [17] states that these criteria hold 
for all but finitely many 7 on the unit circle. 

The start system S{q) = can, for example, be taken 
to be 
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where di is the degree of the i*'* polynomial of P{q) = 0. 
The system (6) is easy to solve and guarantees that the 
total number of start solutions is Y[i=i ^^^^ solutions 
are nonsingular. 

Each homotopy path, starting at a solution of S{q) — 
at t = 1, is tracked to t = using a path track- 
ing algorithm, e.g., Euler predictor and Newton correc- 
tor methods. There are a number of freeware packages 
well-equipped with path trackers such as PHCpack [19], 
HOM4PS2 [20], and Bertini [21]. We used the latter one 
to get the results in this paper. Tracking the solutions 
to t = 0, the set of endpoints of these homotopy paths is 
the set of all solutions to P{q) = 0. Since each homotopy 
path can be tracked independently, NPHC is inherently 
parallelizable. 

The set of real solutions can be obtained from the set 
of complex solutions by considering the imaginary part 
of the solutions (typically, up to a numerical tolerance). 
We remark that the approach of [22] implemented in al- 
phaCertificd [23] can be used to certify the reality or 
non-reality of a nonsingular solution given a numerical 



approximation of the solution. The ability to compute 
all complex solutions, and thus all real solutions, dis- 
tinguishes the NPHC method from most other methods. 
Due to the power of the NPHC method, it has recently 
found several applications in theoretical physics [24]. 

To find the stationary points of the nearest neighbor 
(f)^ model, we need to solve its stationary equations, i.e.. 



= 
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with q^ = {ql,...,q%) G C^. Since (7) is a system of 
coupled third-order polynomial equations, the clas- 
sical Bezout bound is 3^. For this particular system, 
we know that the number of solutions is exactly 3^ 
(counting multiplicity) for any parameters J and /x^ with 
A 7^ 0. This follows since the system consisting of all the 
terms of degree three is a decoupled system of monomi- 
als. That is, there is only one term of degree three for 
the ith polynomial in (7) which depends only upon gf, 
namely the monomial ^{qf)^- This implies (7) has no 
solutions "at infinity" so that the classical Bezout bound 
must be sharp (counting multiplicity). Thus, we have a 
solid check on our claim to find all solutions using ho- 
motopy continuation. However, the problem is that 3^ 
grows rapidly as increases and, due to current compu- 
tational limitations, we are restricted to only small size 
lattices such as 3 x 3 and 4x4. 

For the 3x3 lattice, it took an average of roughly a 
minute to compute the 3^ solutions (counting multiplic- 
ity) for a given value of J using Bertini running on a 2.4 
GHz Opteron 250 processor with 64-bit Linux. For the 
4x4 lattice, it took an average of roughly 8.5 hours to 
compute the 3"'^^ solutions (counting multiplicity) for a 
given value of J using Bertini running on a cluster con- 
sisting of 12 nodes, each containing two 2.33 GHz quad- 
core Xeon 5410 processors running 64-bit Linux. 



IV. PROPERTIES OF STATIONARY POINTS 

Using the NPHC method as explained in the previous 
section, we can obtain all complex stationary points of V. 
In the context of energy landscape methods, one is usu- 
ally interested in the real solutions only, i.e., solutions 
of (7) with q^ G R^. In the next few subsections, we 
report on the properties of these real stationary points: 
In Sec. IV A the number of real stationary points is ana- 
lyzed and the existence of singular solutions is discussed. 
In Sec. IV B we study the potential energies V{q^) of the 
real q'^, and in Sec. IV C their Hessian determinants. In 
Sec. IV D the Euler characteristic of certain submanifolds 
in configuration space, computed from the indices of the 
real stationary points, is investigated. Since, as men- 
tioned in the Introduction and discussed in a Letter [9], 
we found that the real stationary points are not related 
to the phase transition of the model (at least not in the 
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FIG. 1. (Color online) The number of real stationary points 
of V for 3 X 3 (left) and 4x4 (right) lattices, plotted loga- 
rithmically as a function of the coupling J. 



direct way predicted by the theorem in [4] ) , we extended 
our analysis to include complex stationary points. The 
results of this analysis are reported in Sec. IV E. 



A. Real stationary points 

For J = 0, i.e., in the absence of coupling, the sta- 
tionary points of the potential V in (1) can be calcu- 
lated analytically, obtaining 3 ^ disti nct solutions = 
{ql...,q%) with g| € {0 , ± y^eJFJX} . Since A, ^ > 0, 
these stationary points are all real. Upon increasing the 
coupling constant J, real stationary points start to bifur- 
cate into complex ones, and the number of real stationary 
points decreases gradually from 3^ for J = to only 3 
stationary points for some sufficiently large J. This be- 
havior is illustrated for 3 x 3 and 4x4 lattices in Fig. 
1. The three stationary points that persist at large J 
are the two global minima q^ = {ql , . . . , q^j^) where aU 
g| = a/B^^/A, respectively —^/Ofj^JX, and a stationary 
point of index 1 where all qj = 0. 

The value of J at which the number of real solutions 
drops to 3 can be computed semi-analytically by com- 
puting with Mathematica the index of the stationary 
point g** = (0, . . . , 0) as a function of J and then search 
for the value of J at which the index drops to 1. As re- 
ported in Fig. 2, this value is found to be iV-dependent 
and a parabola provides an excellent fit to the numerical 
data of J{N). 

We have also investigated the values of J for which the 
system has at least one real singular solution, i.e., bifur- 
cation points of the parametric systems, using NPHC. 
At these solutions the potential has degenerate critical 
points, a feature that does not make V qualified to di- 
rectly apply Morse theory as described in Section IV D. 
There are two approaches that we used to compute where 
the bifurcations in a one-parameter system occur, which 
we describe in the context of computing where the first 
bifurcation occurs. In the first approach, we use the ba- 
sic philosophy of the NPHC method with a slight change 
that we treat J itself as a continuation parameter, i.e., 
we start with the known solutions at J = and simply 
track the solutions as J increases to determine the small- 
est value of J > where solutions coalesce. This yielded 
the values of J w 0.12907 and J w 0.12894 for the 3 x 3 



and 4x4 lattice, respectively. 

In the second approach, we use the fact that the Hes- 
sian determinant, det7iv{q, J), where. 



(8) 



is zero at the singular solutions. We add this equation, 
det Hvi^q, J) = 0, as an additional equation in the system 
of stationary equations leaving J unfixed so that it can 
be treated as a variable. We then use Bertini to compute 
the set S of values of J where this combined system has 
a solution. Since all of the solutions at J = are non- 
singular, this must hold in a neighborhood containing 0. 
Based on the construction of S [25], we know there must 
be a nonzero univariate polynomial s{x) such that S is 
exactly the set of roots of s. 

The coefficients of the polynomial s depend upon A 
and p?' . If A and p?' are rational numbers, then s has 
rational coefficients meaning that 5 is a finite subset of 
the set of algebraic numbers, a countable subset of C. 
For example, with A = 3/5 and /i^ = 2, we know that 
the set V of complex stationary points must contain 3^ 
distinct points when J is a transcendental number, e.g., 

J = TT. 

For the 3x3 lattice, Bertini found that S consists of 
1357 complex numbers, of which 297 are real and 178 are 
positive. The smallest positive value using this approach 
is also J w 0.12907. This computation also yields that, 
for J > 11.00169, all stationary points must be nonsin- 
gular. Performing this same computation using the 4x4 
lattice is currently beyond the available computational 
resources. 



B. Stationary values 

In the Introduction, we briefiy reviewed the research 
efforts aiming at establishing a relation between phase 
transitions and stationary points of the potential energy 
function V . These efforts all have in common that they 




FIG. 2. (Color online) The value of J at which, for a given 
linear system size L, the number of real stationary points 
of V drops to 3. The dots are data points computed with 
Mathematica, the hue is the parabola 0.0507366L^ fitted to 
the data. 
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FIG. 3. (Color online) The scaled Hessian determinant D plotted vs. the stationary values for all real stationary points of 
a 4 X 4 lattice with couplings J = 0.1, 0.2, 0.3, and 0.5 (from left to right). The distribution of potential energies also illustrates 
that t;'' < for all g'', as discussed in Sec. IV B. 




focus on a conjectured relation between the occurrence 
of a phase transition at some critical potential energy v^. 
and the properties of stationary points with station- 
ary values — V{q^)/N coinciding with Vc [26]. From 
the stationary points obtained by means of the numerical 
homotopy continuation method for lattice sizes 3x3 and 
4 X 4, it is straightforward to compute, via (1), the sta- 
tionary values . For arbitrary couplings J, we found 
that < for all stationary points q^. An analytical 
calculation, reported in [9], has confirmed this observa- 
tion and extended it to lattices of arbitrary sizes. As 
explained in this same reference, it is this upper bound 
on which disproves the theorem by Franzosi and Pet- 
tini [4], as it cannot be reconciled with the fact that the 
critical energy Vc of the phase transition becomes positive 
for couplings J > 0.8. 

C. Hessian determinant 

Once a relation between stationary points of the poten- 
tial energy landscape and the occurrence of phase transi- 
tions had been conjectured in the 1990s, it immediately 
became clear that not all stationary points induce phase 
transitions. Therefore an obvious question to ask was: Is 
there a certain property of a stationary point that ren- 
ders it capable of inducing a phase transition? Some 
years later it was noticed that the Hessian determinant 
Hv of the potential energy function V, evaluated at the 
stationary points, is crucial for discriminating whether or 
not a stationary point can induce a phase transition in 
the thermodynamic limit [27]. For some models, even in 
the absence of an exact solution, this insight facilitated 
the exact analytic computation of transition energies [28] . 
We refrain here from stating the precise criterion, noting 
only that stationary points with a Hessian determinant 
approaching zero in the thermodynamic limit play an im- 
portant role. 

We evaluated the determinant of the Hesse matrix (8) 
at all of the real stationary points q^ of V obtained by 
the homotopy continuation method. In Fig. 3, we show 
the rescalcd Hessian determinant 

D = \dctHv{q')\^^'^, (9) 
plotted versus the stationary values — V{q^)/N for 



all real stationary points of 4 x 4 lattices and various 
couplings J. From these plots one can immediately verify 
that < for all real stationary points and arbitrary 
coupling J, as discussed in Sec. IV B. We will come back 
to these plots in Sec. V when we will use them to compare 
the homotopy continuation data to those obtained by 
means of the Newton-Raphson method. 

D. Euler characteristic 

In the Introduction, and also at the beginning of Sec. 
IV, we referred to the work of Franzosi and Pettini [4] 
or to related publications as dealing with the relation of 
stationary points of the potential energy function V to 
thermodynamic phase transitions. Although this is cor- 
rect as regards content, it is not obvious at first glance, 
as these results were originally phrased in terms of topol- 
ogy changes of certain submanifolds Af„ in configuration 
space r, 

A/„ = {g e r I V{q) < Nv] . (10) 

Upon variation of the parameter v, the topology of the 
submanifolds My may change at some value vt, in the 
sense that My is not homeomorphic to My, for v < vt 
and w > Vt- The occurrence of phase transitions at some 
critical potential energy Vc was then conjectured to be 
related to the presence of topology changes with energies 
Vt in an open neighborhood of Vc- Via Morse theory, 
such topology changes can be related to the presence of 
stationary points of V with stationary values v'^ = Vt (see 
[3] for an elementary introduction). 

In the context of configuration space topology, the Eu- 
ler characteristic xi^^v) of the manifolds Aly has been 
used in several publications as a way of characterizing 
the changes of topology [8, 29]. The Euler characteris- 
tic X is a topological invariant, i. e., different values of x 
for manifolds My and My, imply that My and My, are 
not homeomorphic. Hence monitoring the Euler char- 
acteristic of the family {My}^^^ of configuration space 
subsets under variation of the parameter v, we may get 
an impression of the way the topology of the My changes. 
Plotting the related quantity 

^(^)= J™ i^ln\x{My)\ (11) 
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FIG. 4. (Color online) Graphs of the logarithm of the Euler characteristic, \n\x{Mv)\/N , for a 4 x 4 lattice and coupling 
strengths J = 0.1, 0.15, 0.3, and 0.45 (from left to right). 



as a function of the potential energy v, a kink in a was 
observed precisely at the critical energy Vc of the phase 
transition for several models studied [8, 28, 29]. 

Knowing all stationary points of V with stationary val- 
ues up to a given value v, the Euler characteristic of 
My can be calculated by means of the formula 

N 

x{My) = ^(-l)>,(z;), (12) 

4=0 

where the Morse numbers ^ii{v) are defined in this con- 
text as the number of stationary points of V with index 
i and stationary value < v. The index i is defined as 
the number of negative eigenvalues of the Hessian matrix 
Tiviq^), which is assumed to have only nonzero eigenval- 
ues. As we noted earlier, for finitely many values of J the 
corresponding systems of equations indeed possess singu- 
lar solutions. Using the NPHC method, we know which 
of the values of J possess at least one singular solution 
and in this section, we avoid such values of J. 

We have computed the Euler characteristic x(-^ti) 
from the real stationary points of V as obtained by 
the homotopy continuation method, and the results are 
plotted as a function of v and for various values of J in 
Fig. 4. Since the energy levels are very closely spaced, it 
is difficult to distinguish one from another. Here, we use 
the tolerance 10~^, i.e., if \v2 — f2| > 10~®, then vi and 
1)2 are distinct energy levels. It is anything but a sur- 
prise that no kink or other signature is visible in x(-M„) 
at w = Vc'. As was discussed in Sec. IV B, the stationary 
values are nonpositive, and the Euler characteristic 
x{Mv) is therefore constant for w > 0. In fact, the v vs 

In |x(-M„)| plots do not even exhibit smooth curve like 
behavior as claimed in [8] . Since the critical potential en- 
ergy Vc can be positive for sufficiently large J, it is clear 
that in general it cannot be accompanied by a signature 
in x(M„). 

We can, however, use the computation of the Euler 
characteristic as a consistency check: For potential en- 
ergies V > 0, the manifold My is homeomorphic to an 
-/V-dimensional ball, and the Euler characteristic is there- 
fore known to be x(-^u) = 1 for > 0. Computing 
the alternating sum (12) with all the stationary points 
and their indices as an input, we find that at w = 0, 
x{My) = 1. Since there is no stationary point for v > 0, 
x{My) = 1 for all w > 0. We have confirmed this result 



for all the values of J without singular solutions used in 
this paper. 

E. Complex stationary points 

In Sec. IV B, we discussed the fact that, for arbitrary 
coupling J, the stationary values v'^ are never positive, 
while the critical energy Vc of the phase transition of 
the nearest-neighbor model becomes positive for suf- 
ficiently large J. A direct relation between phase tran- 
sitions and stationary points of V (in the spirit of the 
one in [4]) is hence ruled out, but one might wonder if 
a modification of the conjectured relation might be more 
successful. 

One possible and rather straightforward generalization 
of this conjecture is obtained by considering not only real 
stationary points, but also complex ones. The reasoning 
behind this generalization is that the presence of com- 
plex stationary points whose imaginary parts go to zero 
with increasing system size N should have the same (or 
at least a similar) effect on the thermodynamic proper- 
ties of the system as their real counterparts. To test this 
idea, we have used the (in general complex) stationary 
points q^ obtained by means of the homotopy continu- 
ation method and plotted in Fig. 5 real and imaginary 
parts of the (complex) potential V{q^) for various values 
of the coupling J. 

At first sight the results are encouraging, as they show 
that, for sufficiently large J, there exist complex q^ with 
positive real stationary values V{q^). Moreover, for the 
couplings J we studied, the maximal real stationary value 
is larger than the critical potential energy of the phase 
transition. Unfortunately, from the data we have there is 
not much more we can say, and it would be unreasonable 
to conjecture a relation of the above mentioned kind on 
the basis of our results. 



V. NEWTON-RAPHSON METHOD 

The Newton-Raphson method is a powerful and fre- 
quently used iterative algorithm for approximating the 
roots of a function (see Sec. 9.7 of [30]). In the context 
of energy landscapes, the stationary points of V are de- 
termined by the system of N equations (7), so the prob- 
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FIG. 5. (Color online) Imaginary vs. real parts of the complex potential energies V{q^)/N for all complex solutions g'' of a 4 x 4 
lattice with coupling J = 0.1, 0.2, 0.3, and 0.5 (from left to right). For small couplings J < 0.2 the real part is nonpositive for 
all q^, whereas for larger couplings some of the stationary values move into the right halfplane. 



lem is equivalent to finding the roots of the vector-valued 
function on the left-hand side of (7). 

From a given initial point in phase space, the Newton- 
Raphson method iteratively finds approximations to a 
stationary point. If the function has more than one sta- 
tionary point, it will depend on the initial value of the 
iteration which of the stationary points is found. For 
the potential energy function (1) of the two-dimensional 
nearest-neighbor (p'^ model, we have seen in Sec. Ill that, 
at least for small coupling J, the number of stationary 
points is exponentially large in the number N of lattice 
sites. The result of the Newton-Raphson computation 
will therefore crucially depend on the set of initial points 
chosen for the iterations. First, the initial points have to 
differ sufficiently from each other, in order to make sure 
that different stationary points are found in the various 
iteration runs. Second, properties of the initial points 
will have an influence on the properties of the stationary 
points found, as the outcome of a Newton-Raphson run 
typically yields a stationary point that is in some sense 
close to the initial point. 

For a given coupling J and lattice sizes up to -/V = 
32 X 32, we generated sets of 10^ initial points by means 
of a standard Metropolis Monte Carlo dynamics in con- 
figuration space [31]. The temperature T in the canon- 
ical acceptance rate of the Monte Carlo algorithm was 
set to T = 100, and we will comment on this choice of 
T towards the end of this section. Starting from each 
of the thus generated initial points, the routine newt 
from [30], a globally convergent version of the Newton- 
Raphson method, was used to compute stationary points 
of V. Like in the homotopy continuation computations, 
all stationary points were found to have nonpositive 
potential energies < 0, and the number of stationary 
points was found to decrease dramatically with increas- 
ing coupling J. 

For smaller couplings (J ~ 0.1 and J — 0.2) where 
the number of stationary points is large, we have plotted 
the results of the Newton-Raphson calculations in Fig. 
6. Like for the results from numerical continuation in 
Sec. Ill, we have plotted the scaled Hessian determinant 
Z? at a stationary point versus its stationary value v^. 
For the smaller system sizes N = L x L with L = 3 
and L = 4, the shapes of the clouds of points shown in 
Fig. 6 resemble the ones produced from the complete set 



of stationary points in Fig. 3. For larger system sizes 
L ~ 6, 8, 16, the cloud of points becomes more and more 
focused, being localized in that region of the {v, D) plane 
where the concentration of stationary points is largest. 

We have seen that, in contrast to the homotopy contin- 
uation method where only small system sizes L = 3 and 
L = 4 were accessible, the Newton-Raphson method can 
be applied to much larger sizes up to -L = 32 (and even 
larger with more numerical effort). However, for small 
couplings J and the larger L considered, the number of 
real stationary points of V is expected to be of the order 
of 3^, and it is evident that we can not compute more 
than a small fraction of them. 

This is reminiscent of the situation one encounters in 
Monte Carlo simulations where only a tiny subset of a 
tremendously large configuration space can be sampled. 
In the Monte Carlo context, the problem can be over- 
come (or at least significantly abated) by the technique of 
importance sampling [31]. We have tried a very straight- 
forward (and possibly naive) adaptation of this idea to 
the Newton-Raphson computation of stationary points, 
simply by adjusting the parameter T of the Metropolis 
importance sampling algorithm which was used for gen- 
erating the initial points of the Newton-Raphson search. 
Somewhat disappointingly, the shape of the cloud of 
points in Fig. 6 turned out to be entirely insensitive to 
changes in T. Using for example a small value of T, we 
would have expected to end up with stationary points of 
lower potential energy on average, but surprisingly this 
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FIG. 6. (Color online) Numerical results from the Newton- 
Raphson method. For system sizes N = L x L with L = 3, 
4, 6, 8, and 16, the scaled Hessian determinant D is shown 
versus the stationary value v^. Up to 10® different stationary 
points per system size have been computed for J = 0.1 
(left) and J = 0.2 (right). 
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was not the case. 

There are other, more involved ways of how one could 
shift the search of stationary points to higher or lower po- 
tential energies, but we have not yet implemented such 
refinements. One could, for example, use a more ad- 
vanced search routine (like the OPTIM program package 
[32]) which allows one to search for stationary points of 
a given index, i.e., of a given number of negative eigen- 
values of the Hessian at the stationary point. Since the 
index of a stationary point and its potential energy are 
expected to be correlated, such a routine should find sta- 
tionary points of low energy when searching for small 
indices, and vice versa. 

VI. CONCLUSIONS 

Two numerical methods for the computation of sta- 
tionary points of multivariate functions were discussed in 
this article: the numerical polynomial homotopy contin- 
uation method (NPHC) and a globally-convergent vari- 
ant of the Ncwton-Raphson method. We applied both 
methods to the potential energy function V of the two- 
dimensional nearest-neighbor cj)'^ model on L x L square 
lattices. The NPHC method allows one to obtain all sta- 
tionary points of V, but is limited to system sizes up to 
4x4 with the computational resources we had at our 
disposal. With the Newton-Raphson method we have 
computed stationary points for larger lattices of up to 
32 X 32 sites, but only a small subset of all the stationary 
points of such a large system could be obtained. 

The motivation for this type of study originates from a 
number of conjectures relating the stationary points of V 
to the occurrence of phase transitions in the thermody- 
namic limit. These conjectures refer to certain quantities 
which can be computed from the stationary points of V, 
like their potential energies, their Hessian determinants, 
and the Euler characteristic of the underlying potential 
energy manifolds in configuration space. We have cal- 
culated these and a few other quantities from the sta- 
tionary points of the cj)'^ model obtained with NPHC and 
Newton-Raphson, but — contrary to what the conjectures 
suggest — no sign of the phase transition of the model was 
found. This failure and its consequences, including the 
falsification of a theorem allegedly proven in [4] , was dis- 
cussed in a Letter [9]. 

The NPHC results for the nearest-neighbor (f>'^ model 



on a 4 X 4 lattice can be overviewcd as follows: 

1. The number of real stationary points decreases 
from 3^ for J = to only 3 with increasing J 

2. Singular solutions occur only for finitely many val- 
ues of J. 

3. The stationary values are all nonpositive for ar- 
bitrary couplings J. 

4. The Euler characteristic, computed as the alternat- 
ing sum of the Morse numbers, confirms the cor- 
rect and complete computation of all the stationary 
points. 

5. Unlike real stationary points, complex stationary 
points of V can have positive stationary values, but 
we were unable to identify a relation between these 
positive values and the positive phase transition en- 
ergy of the (f>^ model for larger J. 

Since the Newton-Raphson method yields only a sub- 
set of all the stationary points, we compared these results 
for system sizes up to 16 x 16 to those obtained by the 
NPHC method for 4 x 4 lattices. For this comparison 
we chose plots of the rescaled Hessian determinant D as 
defined in (9) vs. the potential energy v. A comparison 
of different lattice sizes is of course problematic, but a 
general trend can be deduced: For system sizes 8x8 
and larger, the number of stationary points becomes in 
general so large that only that region in the (I?, u)-plane 
is explored where the (strongly peaked) density of sta- 
tionary points is the highest. Importance sampling may 
provide a way out of these difficulties, but we have not 
yet implemented such a scheme. 
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